%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 蚁群算法
% 简介：蚁群算法寻找最优路径
% 作者：Zhaojiang
% 日期：2023/10/12
% 企鹅：277746470
% 最短长度：102.6600
% 最短路径：1->23->31->16->19->10->25->22->2->27->18->15->12->11->34->4->13->
%          20->17->24->14->5->21->8->29->7->9->30->28->26->3->6->32->33->1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;clear;clc;
%% 数据加载预览
tsp=TSP(importdata('./city_data.mat'));
figure('Name','蚁群算法');
tsp.draw_path;
disp 当前路径
tsp.print_path;
pause(2);
close all;
figure('Name','蚁群算法','Position',[200 200 1280 480]);
%% 设定算法参数
rng(1234);
ant_num = 100;% 蚂蚁数量
alpha = 2;% 信息素重要程度因子
beta = 2;% 启发函数重要因子
rho = 0.1;% 信息素挥发因子
Q = 1;% 信息素总量
max_iteration = 1000;% 最大迭代次数
%% 算法数据初始化
city_num = tsp.city_num;% 城市数量
ant = zeros(ant_num,city_num);% 生成蚁群
tau = ones(city_num,city_num);% 蚁群信息素初始化
ant_length = tsp.path_length * ones(1,ant_num);
delta_tau = zeros(city_num,city_num);
best_iteration = 0;
best_path = tsp.path;% 最佳路径
best_length = tsp.path_length;% 最佳长度
all_best_length = zeros(1,max_iteration);% 记录最佳函数值变化
for i = 1:ant_num
    ant(i,:) =[1 randperm(city_num-1)+1];% 随机初始化种群
end
%% 开始迭代求解
for iteration=1:max_iteration
    for i = 1:city_num
        ant(i,2) =randi([2 city_num]);% 随机初始化种群
    end
    city_index = 2:city_num;
    for i = 1:ant_num 
        for j = 3:city_num
            tabu = ant(i,2:(j-1));
            allow_city = ~ismember(city_index,tabu);
            allow_city = city_index(allow_city);
            P = zeros(size(allow_city));
            for k = 1:length(P)
                P(k) = tau(tabu(end),allow_city(k))^alpha*...
                    tsp.distence(tabu(end),allow_city(k))^(-beta);
            end
            P = P/sum(P);
            P_sum = cumsum(P);
            index=find(P_sum>rand);
            next_city = allow_city(index(1));
            ant(i,j) = next_city;
        end
        tsp.update_path(ant(i,:));
        ant_length(i) = tsp.path_length;
        if ant_length(i) < best_length
            best_path = ant(i,:);
            best_length = ant_length(i);
            best_iteration = iteration;
        end
    end
    all_best_length(iteration) = best_length;
    for i = 1:ant_num
        for j = 3: (city_num-1)
            delta_tau(ant(i,j),ant(i,j+1)) = ...
                delta_tau(ant(i,j),ant(i,j+1))+Q/ant_length(i);
        end
        delta_tau(ant(i,city_num),ant(i,2)) = ...
            delta_tau(ant(i,city_num),ant(i,2))+Q/ant_length(i);
    end
    tau = (1-rho)*tau + delta_tau;
    if mod(iteration,max_iteration/100)==0 || iteration ==1
        subplot(1,2,1);hold off;
        tsp.update_path(best_path);
        tsp.draw_path;
        disp(['迭代次数：' num2str(iteration) '/' num2str(max_iteration)])
        disp(['最短长度:' num2str(best_length)])
        subplot(1,2,2);hold off;
        plot(all_best_length(1:iteration));
        drawnow;
%         pause(0.1);
%         break;
    end
end
%% 处理结果
clf;
subplot(1,2,1);hold on;
tsp.update_path(best_path);
tsp.draw_path;
% 打印结果
disp(['最优迭代次数：' num2str(best_iteration)])
disp(['最短长度:' num2str(best_length)])
disp 最佳路径:
tsp.print_path();
% 可视化结果
subplot(1,2,2);
plot(all_best_length)
save('./Best_Path_ACO.mat','best_path','-double');